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Diffuse fields in layered media 1 



Abstract 

The spectral decomposition of the elastic wave operator in a layered 
isotropic half-space is derived by means of standard functional analytic 
methods. Particular attention is paid to the coupled P-SV waves. 
The problem is formulated directly in terms of displacements which 
leads to a 2 x 2 Sturm-Liouville system. The resolvent kernel (Green 
function) is expressed in terms of simple plane-wave solutions. Appli- 
cation of Stone's formula leads naturally to eigenfunction expansions 
in terms of generalized eigenvectors with oscillatory behavior at infin- 
ity. The generalized normal mode expansion is employed to define a 
diffuse field as a white noise process in modal space. By means of a 
Wigner transform, we calculate vertical to horizontal kinetic energy 
ratios in layered media, as a function of depth and frequency. Several 
illustrative examples are considered including energy ratios near a free 
surface, in the presence of a soft layer. Numerical comparisons between 
the generalized eigenfunction summation and a classical locked-mode 
approximation demonstrate the validity of the approach. The impact 
of the local velocity structure on the energy partitioning of a diffuse 
field is illustrated. 



PACS numbers: 43.40.Hb,43.20.Bi,43.40.Ph 
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I. INTRODUCTION 



Since its introduction in elastodynamics by Weaver 1 , the diffuse field concept and the 
principle of equipartition of elastic waves have been successfully applied to both field 2 and 
laboratory experiments^. It is an important ingredient of the spectacular reconstruction of 
Green's function from thermal noise 4 or seismic coda waves^. The principle of equipartition 
can also be used to predict the partition of energy in diffuse fields^ 2 -' 6 . Yet, the practical im- 
plementation of the equipartition principle in a configuration as simple as a half-space poses 
some serious technical difficulties 6 . The simplest mathematical formulation of equipartition 
relies on the the existence of a complete normal mode basis. In the past, simple waveg- 
uides or a homogeneous half-space have been considered 6 - 1 ^. In this work, we construct a 
complete set of normal modes of the elastic wave equation from the general spectral the- 
ory of self-adjoint operators in Hilbert space 8,9 . A diffuse field is then represented as a 
white noise distributed over the complete set of normal modes independent of the nature 
of the spectrum (discrete/continuous), or medium (open/closed). The paper is organized 
as follows: in section [Til the spectral problem is introduced and general properties of the 
elastodynamic operator are summarized. In section IHH the key mathematical results are 
presented and illustrated on a simple problem. Sections IIVIM present the central results of 
the paper. The spectral theory of a layered elastic half-space is derived in detail and applied 
to energy partitioning of diffuse fields. Other applications of the theory are briefly discussed 
in conclusion. 



II. PROBLEM STATEMENT 



The elastodynamic equation governing seismic wave propagation can be written as 



df \u) = -L \u) , 



(1) 
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or more explicitly in the the position representation: 

1 



d t Ui{*) 



djCijkidkUiix). 



(2) 



In equation (TO, the minus sign has been introduced in order to deal with a positive definite 
operator. 

In an isotropic medium where properties depend only on the depth coordinate z, the 
elastodynamic equation decouples into independent scalar (SH) and vectorial (P-SV) equa- 
tions. Our goal is to obtain a complete set of normal modes for the latter problem. Following 



Ref 
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we reduce the three-dimensional problem to a coupled system of second order differ- 
ential equations. Translational invariance suggests to look for solutions of the form u(z)e lkx . 
Additionally, cylindrical symmetry enables one to work in a single plane of propagation 
and ignore the third space dimension. Taking all the symmetries into account yields an 
eigenvalue problem for a second-order differential operator Lr^ 

-d z (z\t\u) + kB l d z (z\u) + k 2 C (z\u) = A (z\u) 

(3) 



(z\h\u) 



p{z 



(z\t\u) = Ad z (z\u) + kB (z\u) 

The operator r provides the tractions generated by the displacement field \u) and p denotes 
the density. The matrices A, B and C are defined as 



A = p{z) 



( a\z) ^ 



V 



o 



c 



f (5\z) ^ 



B = p{z) 



( 2p 2 (z) -a 2 (z^ 



(3\z) 







/ 



(4) 



V 



a 2 (z) 



J 



where a and (3 denote the longitudinal and shear wavespeed, respectively. The vector \u) has 
components (iu z , u x ), where u x and u z denote the original displacement field. This change of 
variables can be represented by a unitary transformation U and turns the original problem 
into an equivalent, easier problem. Indeed, the multiplication of the vertical component by 
i makes L a real operator which simplifies the calculations in the layered case. Once the 
eigenvectors of the transformed operator have been obtained, those pertaining to the original 



operator are recovered by application of the inverse operation W . As shown previously^^, 
the elastodynamic operator is positive definite and self-adjoint in the Hilbert space with 
scalar product 

(u\v) = dzp(z)ui(z)*Vi(z). (5) 
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but loosely speaking, the operator 



The precise domain of definition of L is given in Ref . 
acts on functions of finite elastic deformation energy. In equation ([3]) and (jSJ) we have 
introduced the Dirac bra-ket notation for its compactness and convenience. Since it will be 
used throughout the paper, we have summarized the main properties of the Dirac formalism 
in appendix lAl 



III. STONE FORMULA AND A SIMPLE APPLICATION 

A. Mathematical results 

The Stone formula 13 is a powerful functional analytic result for self-adjoint operators 
in a Hilbert space. This result is particularly useful for the mathematical formulation of 
scattering theory£&^. In this paper, we will formally apply this formula to the elastic wave 
equation to obtain generalized eigenfunction expansions. The eigenvectors |e) are solutions of 
the equation L |e) = A |e) but do not belong to the domain of the operator, i.e., they are not 
solutions with finite energy. These eigenvectors are required to construct the complete modal 
solutions of the elastic wave equation^. Our work generalizes previous results obtained in 
the case of a homogeneous half-spaced^. Before stating the Stone formula, we introduce 
the resolvent G(A) of the operator L: 

G = (L - AT)" 1 • (6) 

Because L is self-adjoint, its spectrum is a subset of the real axis and therefore the resolvent 
is defined for ImA ^ 0. In position and polarisation space, the resolvent is represented by 
the Green matrix 

(z\(L-\iy 1 \z') = G lJ (\,z,z'), (7) 
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which possesses the Hermitian symmetry 

G ij (\,z,z , ) = G ji (\* 1 z' 1 z)*. (8) 

We first recall the fundamental spectral theorem and the Stone formula in an abstract 
setting. To every self-adjoint operator L one can associate a projection operator valued 

function (measure) Pa having the following properties 13 

+00 

f(L) = J /(A)dP A , (9) 

—00 

The spectral theorem guarantees that the spectral family is orthogonal, diagonalizes 
the operator L and provides the completeness relation. A practical method to construct 
the spectral projector P^ is provided by the Stone formula, which relates functions of a 
self-adjoint operator to the discontinuity of its resolvent across the real axis^ 3 -: 

/ (L) = JL hm J [G (A + ie) - G (A - ie)} f(X)d\. (10) 

—00 

For a positive operator, the integral can be taken from to +00 only because the spectrum 
is a subset of the positive real axis. As a particular case of equation (JTUl) . one can formally 
obtain the completeness relation: 

+00 

I = -^-lim / [G (A + ie) - G (A - ie)] dA. (11) 
2m e->o+ J 

—00 

For an operator with a mixed spectrum, with both discrete and continuous parts, equation 
([TT]) will take the following form: 

1= /dA^|e m (A))(e OT (A)|+^|e n )(e n |, (12) 

m n 

where |e m (A)) denotes generalized eigenf unctions of the continuum -we assume that no 
singular continuous spectrum exists- that are normalized in the following sense: 

(e m (A)|e m ,(A')) = 5 mm ^(A-A / ). (13) 

In equation (fT2|) . we have introduced a discrete index in the continuous part of the spec- 
trum, which corresponds to the possibility that the space of generalized eigenfunction be 
multidimensional. 



B. Generalized eigenvectors in homogeneous full-space 



We consider the system ([3]) with constant parameters aoo, and p^. First we construct 
the resolvent L — AI outside of the real axis. In order to uniquely define the square root of 
a complex number z we adopt the following convention: y/z = \z\ 1 ' 2 e 10 /' 2 , with 9 = arg(z) e 
[0, 27r) . This definition guarantees that lm.y/z > when z lies outside the positive real axis. 
The 2x2 system of coupled second-order equation ([3]) has four linearly independent solutions 
\ u +oo(ty)i | M -oo(^))' | u +oo(^))' \ u -oo(^))i wnose properties are summarized below: 



(^k±oo(A)) 



e a {\) 




(14) 



D ±iq (\)z 



In equation (fT4l) . the prefactors e Qj/ g ensure that the plane waves are energy normalized 10 : 

1 



2p oV / Ag a ,/3 

and the vertical wavenumbers q a ^ are defined as 



(15) 



— - k 2 

PI 



(17) 



'A. = \/ — - } ' 2 - %=\\—- k2 - ( 
The following symmetry relations will prove useful: 

(z\u P +OQ (\l) = - (z\u p +oc (\)) 

(z\u P 00 (X*)) = -(z\u p 00 (X)r 
with completely analogous relations involving |w+oo(A)), | m5 oo(A))- When ImA = 0, A > 
a^/c 2 , assuming a time dependence of the form e~ 4v/ ^', (u^, u^oo)-^?^, u? M ) behave like 
(outgoing, incoming) plane P-S 1 waves at +oo, respectively. Using this set of solutions, we 
can construct the Green matrix for ImA 7^ 0, which obeys the following equation: 

- p^d z r[G)ij(z, z) + K ik (z)G kj (X, z, z) - XG^z, z') = p^5ij5(z - z'), (18) 

where we have introduced the differential operator K(z) = k 2 C + kB t d z . The Green matrix 
should possess the following properties: 



1. For z > z' , each column of G must be a linear combination of solutions with finite 
energy at +00 , |«Joo)>- 

2. For z < z' , each column of G must be a linear combination of solutions with finite 
energy at -00 fli^J , ■ 

3. G obeys the symmetry relation (jSj) 

4. G is continuous on the diagonal z = z'. 

5. The traction matrix of G has a jump discontinuity on the diagonal: T[Gr]y|*v_ = — 5^ 

Properties (l)-(5) are straightforward generalizations to systems of 2 nd order differential 
equations of the standard Sturm-Liouville theory^. Properties (l)-(3) indicate that G(A) 
should have the following form 

^i(A)h P 00 (A))« 00 (Ar|+ir 2 (A)| M s 00 (A))<^ 00 (Ar| z<z' 
G(A) - ^ . (19) 

^i(A) |<oo(A)> ^ P oo(A)*| + K 2 (X) \u s +00 (X)) («f 00(A)* I x > z> 

In equation ( [191) we have separated the Green matrix into a P and an S wave part. 
Application of the continuity condition of the Green matrix on the diagonal (4) yields 
-Ki(A) = K2 (A). Next we apply the condition (5) for the discontinuity of the traction 
on the diagonal to obtain: 

K^X) = K 2 {\) = (20) 

The next step is to study the discontinuity of the resolvent across the real axis. We 
will denote by and 1^+^) the limiting values of the vectors ^^(A)) as A approaches 

the real axis from above and below, respectively, with analogous definitions for \u F J^ Q ) and 
|w-f^o), • • ■ • The jump of the resolvent G + — G~ across the real axis depends on the relations 
among these functions. For z < z', one obtains 

g + - g- = -j= (10 <«+n + \u p sj + > («+ y\ + 10 («-)*!) . 

V A 

(21) 



Three cases have to be distinguished: 

(1) A > a^k 2 . In this case q a $ and e a ^ are real and all the waves are propagating. Applying 
the following symmetry relations, 

K«Xr>=io io=-i(«ar> 

(22) 

K«2or>=i«2o> = - \« + ooT) , 

one can rewrite the resolvent discontinuity fl2Tl) as a sum of projection operators on a gen- 
eralized eigenfunction basis: 

G + - G = -JL= (| M p +> | + I + |« 5 + > I + |<+ ) «+ I) . (23) 



From equation ( l23i) . one concludes that there are four linearly independent generalized 
eigenfunctions with eigenvalue A > a^k 2 . For instance, the first term on the left-hand side 
of equation ( 1231) is an orthogonal projector on the space of upgoing P waves. 
(2 ) a^k 2 > A > P^k 2 . q a is pure imaginary and q@ is real which corresponds to evanescent 
P waves and propagating S waves. Upon application of the following symmetry relations: 



I («-«,)) = ~i Koo) = -* \ U +to) 

K«S,r> = i«s,> = -i(«£,r> 



(24) 



equation ( 1211 simplifies to 



G + " G = -= ) (u'U + \u s + to) (u S + to\) • (25) 
V A 



There are two linearly independent S eigenfunctions. 

(3) A < P^k 2 . There are no singularities and the resolvent is continuous: 

G+ - G = 0. (26) 

Because of the symmetry of the Green's function (151) , the same results are obtained for 
z > z'. The completeness relation now follows from the Stone formula ( TTTI) . The continuous 
parameter A is analogous to a squared eigenfrequency. Introducing A = u 2 , and generalized 



eigenvectors ^^(o;)), ■ • ■ the resolution of the identity writes as a sum of integrals over 
frequencies: 



+00 



+00 



(27) 



du; {\e s +00 (u)) (e s +00 (u;)\ + |e^M> (e^Hl) 



/3oofc 



The eigenvectors ^^(a;)), • ■ ■ can be deduced from the basis vectors |m+oo(A)) and form 
an orthonormal set. Let us for instance consider the eigenvector le+^u;)): 



(*|eJooM) 



00 



\ k I 



D +iq a z 



(e P +00 {u)\e p +00 {u'))=5{u-J) 1 (28) 



1 „ikx 



where q a = \/ ^2 k 2 . Multiplying the wavefunctions (zle+^u)) by the phase term -j^e 



we can form compound eigenvectors |^ P ' ) that can be used to expand 2-D in-plane vector 
fields. The expansion takes a particularly simple form if instead of using a frequency inte- 
gral, the new variables p z = */ — k 2 and p z = \l — k 2 are introduced in the P and S 



integrals of equation (1271) . respectively. In order to obtain eigenfunctions of the elastody- 
namic operator in its standard form, we apply the unitary transformation U* to obtain a 
fairly simple completeness relation: 



I 



dp x dp z (\?P P (PxiPz)) {V(Px,Pz)\ + \4> S (Px,Pz)) (lp S (Px,Pz)\) 



(29) 



where \4> p (p x ,Pz)) an d \i> s (Px,Pz)) are simple plane P and SV waves: 

e ip x x+ip z z 



(r\ip (p x ,p z )) =p- 



(r\?p S (p x ,Pz)) =P A 



Jp x x+ip z z 



(30) 



where p and p 1 - denote unit vectors parallel and perpendicular to the wavevector p, respec- 
tively: 



P 



1 



ip 2 x + p 2 



H 

\Pz) 



1 



p 



Vpx+pI 



\Px J 



(31) 
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\ip p (jpxiVz)) an d \'4 )S {'PxiPz)) are eigenvectors of the elastic wave equation with eigenvalues 
a ooiPx + Py) an d P%c(Px + respectively. They form a complete set of normal modes that 
decompose 2-D vector fields into longitudinal and transverse parts. 



IV. GENERALIZED EIGENFUNCTIONS IN STRATIFIED MEDIA 

A. Construction of the resolvent 

We begin with the construction of the resolvent, which parallels the development of the 
previous section. To make comparison easier we denote by p^, and (3^ the density and 
seismic wave speeds in the half-space located at depth z > 0. The medium is assumed to 
be composed of a stack of n layers and is bounded by a free surface at z — —z n , z n > 
0. We wish to find the eigenfunction expansion associated with the eigenvalue problem 
L \u) = X\u) introduced in equation ([3]) and supplemented with the free surface traction 
condition (z\t\u) = at z = and the continuity of stresses and displacements at each 
interface. To do so, we introduce four linearly independent vector solutions of the elastic 
wave equation (jHJ) denoted by |w+oo(A)), ^^(A)), jw+^A)), |w£oo(A)). In the half-space 
z > 0, their analytical form is given by equation (TH1) . Their analytical dependence in the 
stack of layer is obtained by integrating the equations of motion ([3]) from depth z = + to 
z = —z+ by applying continuity of stresses and displacements at each interface. Note that 



id, 



in general, these solutions do not verify the stress free condition at z = 0. Following Ref. 
we introduce two linearly independent solutions Wq (A)), |mq(A)) of the wave equation that 
verify the stress-free condition at z = 

K(A)> = Koo(A)) + r^(A) KJA)) + H-(A) (^(A)) 
K(A)> = h S oo(A)> + r^(A) \u p +00 (\)) + r"(A) ^(A))' 
where r pp , r ps , r sp and r ss are the reflection coefficients of the stack of layer including the 
free surface. They are uniquely specified by the boundary conditions and the form of the 
solutions in the half-space. For ImA = 0, ReA > a^A; 2 , the solutions ( 1321) behave like 
propagating P and S waves incident from +oo together with their reflections. They are 

11 



solutions of the wave equation with infinite energy that verify the free surface condition. 
For ImA 7^ 0, we remark that because the matrix elements of L are real, the following 
relation holds: L \uq(X)*) = X* \uq(X)) = L \uq(X*)). Using equation (TIT]), which is still 
valid in the stratified case, the following symmetry relations are established: 

|<' 5 (A*)> = - |<< 5 (A)> , (33) 

r ab (X*) = r ab {X)*, (34) 
where a,b = p, s. In addition, using the method of propagation invariants developed in 



RefJl7l. we can prove the following reciprocity relation: 

r P s (X) = r sp (A). (35) 

To construct the Green matrix we again make use of the properties (l)-(3) above. Condition 
(1) is replaced by the requirement that the columns of G be linear combinations of the vectors 
\v,q), \uq) for z < z' . This suggests to try the following form for G(A): 

ir 1 K(A)>« 00 (A)*|+ir 2 |^(A)><^ 00 (Ar| z<z> 
G(A) - { , (36) 

K x \u p +00 {X)) «(A)*| + K 2 \u s +00 {X)) (u s (X)*\ z > z> 

where use has been made of equation ( 1341) . In essence, equation ( 136]) is completely similar to 
equation ( fl~9l) : the Green function is separated into a "P" and an "S"' wave part, although 
the vectors \v.q), \v,q) are neither purely longitudinal nor purely transverse, due to the mode 
conversions at the medium boundaries. We must now determine the values of Ki and K 2 in 
order to fulfill conditions (4)-(5). When the source is located in the half-space z > 0, we can 
make use of the analytical form of the vectors \uq' s ) and \u+^). Following the same steps 
as in the homogeneous case and employing the reciprocity relation (135|) one obtains: 

Kl (X) = K 2 (X) = (37) 

Let us now show that this expression is valid throughout the stack of layers. The key point is 
to verify conditions (4)-(5). Let us first introduce the following 4-dimensional displacement- 

12 



stress vectors : 



and 



\t[u+'^(\, Z)} j 



(38) 



w ' (X,z) 



(39) 



The conditions (4)-(5) will be satisfied provided the following relation holds at arbitrary 
depth z in the medium: 



^= z)w%(\, zf - u£(A, z)w?(X, zf 

+wi(\, z)w$(\, zf - wg(X, z)wi(X, zf) = N, 



(40) 



where we have introduced the notations : 



N 



E 



and E 



1 



v° V 



(41) 



Relation (l4"0"j) is verified in the half-space z > as discussed above. The displacement-stress 
vector satisfies a first-order linear system dw(X, z) / dz = A(X, z)w(X, z), equivalent to the 
second order system (j3J) , which can be integrated with the aid of a propagator matrix: 



w(X, z) = P(X, z, z')w(X, z). 



(42) 



Since the propagator matrix has the symplectic symmetry^ - 



P(X,z,z')NP(X,z,z'f = N, 



(43) 



relation (HOI can be propagated at any depth in the stratified half-space by left and right 
multiplication of equation (140]) by P and its transpose, respectively. Note that in the case 
of a stratified medium, the symmetry relation (1431) can be verified by direct but tedious 
calculations. 
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B. Singularities of the resolvent 



The next step is to reexamine the discontinuity of the resolvent (1361) across the real 
axis. This task is significantly facilitated by the unitary and symmetry properties of the 



reflection coefficients. As shown in RefJl7l. the energy normalization of the eigenvectors 
is of fundamental importance in this respect. Other normalizations are of course possible 
but they lead to much more awkward calculations. We will denote by \u(X + ))-(\u(X~))) 
the limiting value of the ket \u(X)) as A approaches the real axis from above- (below). The 
interrelations among the kets across the real axis determine the form of the 

resolvent discontinuity. These relations, in turn, depend essentially on the analytical forms 
of these kets, i.e., on whether the P and S waves are propagating or evanescent in the half- 
space z > 0. As in the homogenenous case, this leads to a separation of the real axis into 
three distinct domains, as shown in Figure [U 

(1) For propagating P and S waves: A > a 2 k 2 , the following relations apply: 



u 



±00 



u 



±00 



(A+)> = Koc(A + r) hLo(A + )> = KJA+)*) 

(a-)) = - Koo(A + r> hLo(A-)) = - Koo(A + r> 



(44) 



These relations, together with the analytical properties of the reflection coefficients, enable 
us to write the kets \uq ,s (X~)) in terms of the kets |w^(A + )) only: 



u 



p (A-)> = - |< 00 (A + )) - ?pp KJA+)) - r^\ui 0O {X + )) , 

(45) 

s (A-)> = - \u s +00 (X + )) - KJA+)) - —v \u^(X + )) , 



where the reflection coefficients are defined by their limiting value as A approaches the real 
axis from above. For notational convenience, we have denoted the complex conjugation 
with an overbar. When both P and S waves are propagating in the half-space, the matrix 
of reflection coefficients is unitary^: 



R 



ir-PP r Ps\ 
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, RR j = E. (46) 



For z < z', the discontinuity of the resolvent writes: 

G + " G- = -j= (K(A+)) <«f 00 (A + )| + K(A+)) <^ 0O (A + )| 
- K(A")> «oo(A + )| - K(A-)) <^ 00 (A + )|) 
Using equations f j4"5l) and ( 1461) and after some algebra, one obtains: 

G + - G- = (|^(A + )> <^(A+)| + \u s (X + )) <^(A+)|) . (48) 



The same result applies for z > z' and we have therefore put the discontinuity of the resolvent 
in suitable dyadic form. It appears that for A > a 2 k 2 there are two linearly independent 
eigenvectors. They are energy normalized incident P and S waves incident from +00 together 
with their reflections. 

(2) For ot 2 k 2 > A > (3 2 k 2 , there are propagating S and evanescent P waves in the half-space 
z > and the following relations apply: 



\ U ±oo 



(A + )) = - i ^(A+r) and \u s ±00 (X + )) = |^ oc (A + )*> , (49) 



which yields: 



u. 



'(A-)) = ~i Koo(A + )) - i?* Koo(A + )) - ^ h£oc(A + )> 



?(A-)> = " h+oo(A + )> - |^ 0O (A + )) - |< 00 (A+)> 



(50) 



For propagating S 1 and evanescent P waves in the half-space, the unitary relations take 
a more complicated forrn^: r ss r sp = ir sp . With the aid of equation (1491) . this implies 
r ss |mq(A + )) = — |uq(A - )) . For z < z', using the unitary relations, the resolvent disconti- 
nuity can be put into the following form: 

G + -G- = -[ ? |^(A + )) + K(A-))]« 00 (A + )| 

(51) 

+ \u s (\ + ))[(u s _^\+)\+^(u s +oo (\ + )\] 
Using the symmetry relation r pp — r pp = i r sp , the following relation is established: 

i k P (A + )> + \u^X-)} = -~ p |4(A+)> , (52) 
15 



which in turn implies: 

G + - G" = |«?(A + )> <uf(A+)| . (53) 

This last equation has the desired form and shows that in the range of A considered there 
is only one generalized eigenfunction. In the half-space, it consists of the sum of incident, 
reflected propagating S wave and a reflected evanescent P wave. Note that for some A 
in the range }P^ k 2 ,al k 2 [, it may happen that the vector (w^-,) satisfies the zero traction 
condition at the free surface. Such a situation can be considered as accidental, since a 
slight modification of the thickness or velocity in one the layers is likely to make the mode 
disappear (Y. Colin de Verdiere, personnal communication) . Such possible complications 
will be neglected. 

(3) For j3 2 k 2 > A, both P and S waves are evanescent in the half-space z > 0. In this case 
the reflection matrix R is simultaneously symmetric and hermitian, so that all its coefficients 
are real. Further, using the relations: 

l«L,(A-)> = - i l«Lo(A + )> and |4oo(A-)> = - i kLo(A + )> , (54) 
one establishes: 

|t£' ff (A-)> = -i |«? 5 (A + )> , (55) 

which seems to implies that the resolvent discontinuity vanishes identically. This is not 
true as there can exist values of A for which the reflection coefficients have poles. For 
these particular points of the spectrum, the kets \uq' s (\ + )) can be written solely as linear 
combinations of the kets |u+'oo(A + )}^- This means that there exist eigenvectors with finite 
energy that satisfy the traction free condition, which are the well-known Rayleigh surface 
waves. To each pole, we can associate a one dimensional eigenspace spanned by a normalized 
Rayleigh wave eigenvector. For a given value of the horizontal wavenumber k, there can be 
n distinct solutions denoted by Ao < Ai < • ■ • < A n . The index serves to label the different 
surface waves mode branches. The index refers to the fundamental Rayleigh mode. 
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FIG. 1. Spectral properties of the elastodynamic operator in a layered half-space. A is the 
eigenparameter. 

C. Eigenfunction expansions 

We have now identified all the singularities of the resolvent on the real axis and are in 
position to write down the completeness relation for 2-D in-plane problems using formula 
(II Op . In the frequency domain, introducing the new variable u 2 = A + , one obtains: 

+oo +oo \aaok\ 

1 = J dkJ2\ eR ( UJ n))(e R (co n )\ + J dk J du\e s {u)){e s {u)\ 



oo +oo 

dk 



-oo [A»fc| 
duo (\e s (uj)) (e s (uj)\ + \e p (uj)) (e p (uj)\) 

-oo \ aoo k\ 

where the compound generalized eigenvectors le 5 '^) are defined as: 



(56) 



(r\e s ' p (oj)) =v / 2^ (z\u S > P (lu 2 )) and (r\e R (uj)) 
\/2n 



Akx 



[z\u R (uj n )). 



(57) 



In the first of equations (1571) . the a/2 factor is the result of introducing of the new variable 
uj, whereas the last equation defines properly surface waves with unit-energy normalized 
z-eigenfunction. Although this is implicit in equation (|57|) . the reader should keep in mind 
that all eigenvectors are functions of the horizontal wavenumber k. This applies to the 
eigenfrequencies u n as well as the number of branches. As in the case of the homogeneous 
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space, it is possible to replace the frequency integral by a vertical wavenumber integral in 
equation f[5"6T) to obtain: 



+00 +00 

I 



J dp x ^2 \ip R {uj n )) {ip R {u: n )\ + J dp x J dp z \ij s (p x ,p z )) (?p s (p x ,p z )\ 



hoo 



-00 -00 (5g) 



+ dp x dp z \tfj p (p x ,p z )) (ip P (p xjPz)\ 



-00 —00 



The eigenvectors \ip p,s ' R ) are obtained by application of the transformation U* to the kets 
\e s,p ' R ). The eigenvectors |^ p,s ) correspond to properly normalized incoming P and 5" waves 
exactly as given in equations (l30]) - fl3Tj) together with their reflections from the stack of layers. 
This is a general result in scattering theory: the eigenvectors of the medium+scatterer 
can be obtained by calculating the complete response of the scatterer to the unperturbed 
eigenvectors^ 1 ^. This is basically what the double integrals of equation (158]) say. Note that 
by calculating the scattering of plane waves, one does not obtain directly the surface waves. 
They show up as poles of the reflection coefficient located on the real axis and, in that sense, 
can be compared to the bound states in scattering theory. 



V. APPLICATION TO DIFFUSE FIELDS 

A diffuse field is defined as a random field where all the modes are equally represented. 
More precisely, according to RefQ, it is a narrow-band signal which is a sum of normal 
modes of the system excited randomly at equal energy: 

u(x,f) = X; // ^^A^^pyM^^^e^^'^ 

(59) 



n,R,L cR 



+ J d Pz Y,JJ^^A p > s (p x ,p y ,p z ^ p f( z y^^ 

p z <0 P ' S ujGB 

In equation (|59|) . B denotes the frequency band of the signal and u refers to the complex 
analytic signal whose real part is the measured displacement vector. In a diffuse field, the 
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amplitudes A P,L and A p,s are assumed to be white-noise random processes: 

{A n,L A P, S) =0; 

M L ( Px ,Py)A% L (p' x ,p' y y) =a 2 5 R > L 5 nm 5(p x -p' x )5(Py-P y )i ( 60 ) 

(A P ' S (p x ,p y ,p z )A P ' S (p',,Py^y) =^ PS S(p z -p' z )5( Px -p'MPy-Py)- 

In the case of a semi-infinite medium, the sum over all modes involves a continuous index 
that labels the eigenfunctions of the continuous spectrum and a discrete index that refers 
to the surface wave mode branch. 

Formulas (!59|) - (l60l) are now applied to the calculation of the vertical to horizontal kinetic 
energy ratio of a diffuse field in a stratified half-space. The vertical (or horizontal) kinetic 
energy density at central frequency uq is obtained by means of the Wigner distribution of 
the complex analytic wavefield: 

E z {t, r, x) = ±p(z) (d t u z (t + r/2, x)d t u z (t - r/2, x)*) . (61) 

In equation (IBTj) . the brackets denote an ensemble average. Inserting the spectral decompo- 
sition ( 159|) . applying the equipartition principle ( 1601 . and taking a Fourier transform with 
respect to the time variable r, the local vertical kinetic energy density of a diffuse field at 
circular frequency ujq is obtained: 



E z (u ,z) = 




In equation (1B"21) . the delta functions represent the density of states of the surface and body 
waves, which are closely related to the imaginary part of the Green's function^. The link 
between the Wigner distribution of the wavefield and the Green's function is the theoretical 
basis of the Green's function reconstruction in diffuse fields^ 1 ^ 1 ^. Introducing cylindrical 
and spherical coordinates in the double and triple integrals of equation fl62|) . respectively, 
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FIG. 2. Depth dependence of the vertical to horizontal kinetic energy ratio near the free 
surface of a Poisson half-space. Dots: locked-mode approximation. Dashed line: generalized 
eigenf unctions summation. The depth unit is the shear wavelength \ s 



In equation (163|) . we have introduced the following notations: c n and u n are the phase and 
group velocities of the n th surface wave mode, respectively; v p -(v s ) stands for Qoo-(Poo)- The 
calculation of the eigenfunctions and the remaining sum and integral have to be performed 
numerically 

In Figure ([2j), we consider the vicinity of a free surface which has been previously in- 
vestigated by several authors^I The calculations were performed both with formula (|63|) 
and with a locked-mode approximation where the medium is bounded at great depth by a 
rigid boundary where the displacements vanish exactly. In the latter case, the eigenvalue 
problem in the depth coordinate only has a discrete spectrum, which is standard. For more 



one finds: 





(63) 



p,s 
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details on the locked-mode technique and geophysical applications, the reader is referred to 



Ref 



22 



The outcome of the two calculations for a 3-D medium are superposed in Figure [2J 
In the locked-mode technique, the lower boundary is at a depth of 16 shear wavelengths. A 
total of 32 Love and 50 Rayleigh modes were found. The agreement is very satisfactory and 
confirms the validity of our approach. Our work shows that the generalized eigenfunctions 

beit with a continuous index. 



using the eigenfunctions 



of the continuum can be treated like standard normal modes, a 
The results presented in Figure (j2J) were first obtained in Refsj2|j7 
of a thick plate, the so-called Lamb modes. Our calculations illustrate the fact that the 
elastodynamic operator is in the limit point case at +oo according to the classification of 
Ref 



, 23l . This simply means that independent of the self-adjoint boundary condition imposed 
at the lower boundary -Neumann, Dirichlet or mixed-, the eigenfunctions converge to a com- 
mon limit as the depth of that boundary tends to oo. A few wavelengths away from the 
free surface, the kinetic energy ratio oscillates around 0.5, which is the expected value in a 
homogeneous half-space. 

We now investigate the case of a soft layer overlying a homogeneous space. In the layer, 
the P and S wave velocity is one third smaller than in the half-space and the density is 
reduced by a factor 2. The vertical to horizontal kinetic energy ratio is plotted as a function 
of depth in Figure [31 The unit depth is the shear wavelength inside the layer and the 
interface between the layer and the half-space lies at about 0.2 unit depth. In the locked- 
mode approximation, two situations were considered where the rigid boundary is located at 
5.25 and 42 wavelengths below the bottom of the layer. In the former case, 7 Love and 11 
Rayleigh modes were found, while in the latter case we identified a total 56 Love and 88 
Rayleigh modes. In the generalized eigenfunction expansion method, only the fundamental 
Love and Rayleigh modes are present. Thus, in the locked-mode method, the higher Love 
and Rayleigh modes serve as an approximation for the propagating P and S waves incident 
from below the layer. As illustrated in Figure [3j the agreement between the two methods is 
extremely good and as expected, improves as the depth of the rigid boundary increases. In 
the vicinity of the layer, the vertical to horizontal kinetic energy ratio shows rapid variations 
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FIG. 3. Depth dependence of the vertical to horizontal kinetic energy ratio in the presence 
of a soft layer. Dots: locked-mode approximation. Dashed line: generalized eigenfunctions 
summation. The depth unit is the shear wavelength X s . In the locked mode approximation, 
the lower boundary is located at a depth of 5.25A S (Top) and 42A S (Bottom), respectively. 
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and, at greater depth, converges towards the ratio of a homogeneous space. 

In Figure HI we investigate the frequency dependence of the V 2 / H 2 kinetic energy ratio 
at the surface of a solid with a soft layer at the surface, where the shear velocity and 
density are reduced by a factor 2 with respect to the underlying Poisson half-space. The 
ratio of longitudinal to shear wavespeeds in the layer is taken to be 2.5. To facilitate the 
interpretation of the V 2 / H 2 calculations, we show in the Top panel the normalized density 
of states of body, Rayleigh and Love waves as a function of frequency. The frequency unit 
is the fundamental resonance frequency of the layer for vertically propagating shear waves, 
fo — (3/4:h, where (3 is the shear wavespeed and h is the layer thickness. In Figure HJ we 
have also indicated high- and low-frequency asymptotics of the V 2 / H 2 energy ratio. At high 
frequency, we expect the waves to be largely insensitive to the properties of the underlying 
half-space. We therefore calculate an approximate high-frequency V 2 / H 2 ratio by replacing 
the true model with a simple half-space where the P to S velocity ratio equals 2.5. The 
expected value is 0.513, in good agreement with the full calculations. At low-frequency, we 
can neglect the presence of the shallow layer to recover the classical 0.56 ratio at the surface 
of a Poisson solid. The V 2 / H 2 ratio is slightly smaller at high frequencies mainly because 
of the increased amount of horizontally polarized shear waves. To the contrary, the V 2 / H 2 
energy ratio of the Rayleigh wave tends to increase with the P to S velocity ratio. Close to 
the resonance frequency fo of the layer, the vertical to horizontal kinetic energy ratio drops 
dramatically. Careful analysis reveals that this is caused by the increasing contribution of the 
fundamental Love mode to the density of states at the surface. Around twice the resonance 
frequency, the V 2 /H 2 ratio presents a marked overshoot which is due to the appearance 
of Rayleigh waves higher modes (i.e. body waves trapped in the low velocity layer) which 
are preferentially polarized on the vertical axis. At high frequencies, the energy density is 
largely dominated by the Rayleigh and Love waves trapped in the low-velocity layer. The 
fluctuations of the local density of states slowly decrease with increasing frequency. 
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FIG. 4. Properties of a diffuse field at the surface of a half-space with a superficial low- 
velocity, low-density layer. Top: Normalized density of states of Rayleigh (dashed line), 
Love (solid line) and bulk waves (dash-dotted line) as a function of normalized frequency. 
The resonance frequency /o of vertically propagating S waves in the layer is the unit fre- 
quency. Bottom: Frequency dependence of the vertical to horizontal kinetic energy ratio. 
Dashed lines: low- and high-frequency asymptotics. Solid line: generalized eigenfunctions 
summation. 
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VI. CONCLUSION 



In this work, we have shown that the definition of a diffuse field as a white noise in modal 
space can be applied equally well to closed- and open systems^. In the framework of the 
spectral theory presented in this paper, discrete, continuous or mixed spectra can be treated 
on the same footing. Other applications of the generalized normal mode expansion could be 
envisaged, such as the calculation of synthetic seismograms for geophysical applications^. 
The theory could also be used to give more rigorous foundations to empirical civil engineering 
techniques. In particular, the large drop of the vertical to horizontal kinetic energy ratio 
in diffuse fields close to the resonance frequency of a low-velocity layer sounds reminiscent 
of the so-called Nakamura's method used for site effect evaluation with ambient noise^. 
Although extremely popular, the limitations of the technique are still to be understood. 
The diffuse field concept offers a potentially useful tool for this purpose. 
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APPENDIX A: SUMMARY OF DIRAC FORMALISM 

In this appendix, we summarize the bra-ket notations used in this paper. A vector or 
ket of an abstract vector space will be denoted by: |e). The representation of the vector in 
position space is written as: (x\e) = e(x). To each ket of our space of function, we associate a 
bra denoted by (e\ = \ey . In position space, a bra has representation (e\x) = e(x)*, where * 
denotes complex conjugation. Mathematically speaking, a bra would be more appropriately 
defined as a linear functional acting on a space of test functions. However, we will not insist 
on these technicalities. We shall also make use of the conjugated versions of kets such that: 
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(x\e*) = e(x)*. The scalar product between two vectors \e) and |/) is denoted by 

(e|/> = J dxp(x)e i {x)*f i (x), (Al) 

where a summation over the repeated index % is implied. The completeness relation or 
resolution of the identity for a set of eigenvectors \e p ) is written as: 

p 

where I denotes the identity in the abstract vector space, and p denotes a label running over 
the whole set of eigenf unctions. Equation (1A2I) introduces the outer product between a bra 
and a ket. The outer product of two (properly normalized) eigenvectors is an orthogonal 
projector on the subspace generated by |e). In the position representation, equation flA2[) 
reads 

^2(x\e p )(e p \x') =$ij — 7^~ - (A3) 
p 

The appearance of the weighting function p in the denominator is consistent with the def- 
inition of the scalar product flAll) . The matrix elements in position and polarization space 
of a general outer product between a ket |e) and a bra (/| are given by 

(x\e)(f\x') = e l (x)f 1 (xy. (A4) 

The matrix elements of an abstract operator L in position space are given by the kernel of 
an integral operator. If L represents a differential operator, it will be assumed to be diagonal 
in position space. 
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